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ABSTRACT 



Context. The migration of low-mass planets, or type I migration, is driven by the differential Lindblad torque and the corotation torque 
in non-magnetic viscous models of protoplanetary discs. The corotation torque has recently received detailed attention, because of 
its ability to slow down, stall, or reverse type I migration. In laminar viscous disc models, the long-term evolution of the corotation 

\l ■ torque is intimately related to viscous and thermal diffusion processes in the planet's horseshoe region. It is unclear how the corotation 

torque behaves in turbulent discs, and whether its amplitude is correctly predicted by viscous disc models. 

Aims. This paper is aimed at examining the properties of the corotation torque in discs where magnetohydrodynamic (MHD) turbu- 
lence develops as a result of the magnetorotational instability (MRI), considering a weak initial toroidal magnetic field. 
Methods. We present results of 3D MHD simulations carried out with two different codes. Non-ideal MHD effects and the disc's 
vertical stratification are neglected, and locally isothermal disc models are considered. The running time-averaged tidal torque exerted 

ri \ by the disc on a fixed planet is evaluated in three different disc models. 

O „ Results. We first present simulation results with an inner disc cavity (planet trap). As in viscous disc models, the planet is found to 

I . experience a positive running time-averaged torque over several hundred orbits, which highlights the existence of an unsaturated coro- 

O ' tation torque maintained in the long term in MHD turbulent discs. Two disc models with initial power-law density and temperature 

profiles are also adopted, in which the time-averaged tidal torque is found to be in decent agreement with its counterpart in laminar 

C/3 , viscous disc models with similar viscosity alpha parameter at the planet location. Detailed analysis of the averaged torque density 

Cu ■ distributions indicates that the differential Lindblad torque takes very similar values in MHD turbulent and laminar viscous discs, 

and there exists an unsaturated corotation torque in MHD turbulent discs. This analysis also reveals the existence of an additional 
corotation torque in weakly magnetized discs. 

►^ , Conclusions. Our results of 3D MHD simulations demonstrate the existence of horseshoe dynamics and an unsaturated corotation 

I . torque in weakly magnetized discs with fully developed MHD turbulence. 

C^^ . Key words, accretion discs — magnetohydrodynamics (MHD) — turbulence — methods: numerical — planetary systems: planet- 

C^^ ' disc interactions — planetary systems: protoplanetary discs 
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O 1. Introduction I migration (e.g|Ida&Linll2008albl:ISchlaufmanet"ai] 120091: 
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(e.g.. Ilda <3 

iMordasini et al.l l2009h . The latter has been intensively stud- 
From the early stages of their formation, planets tidally interact jg^j ;„ two-dimensional (2D) non-magnetized, laminar viscous 
with their natal protoplanetary disc. The tidal torque exerted by ^isc models. In such models, the torque comprises the differen- 
the disc on a planet drives the planet's radial migration. Several tj^j Undblad torque and the corotation torque. The differential 
migration regimes are usually distinguished, depending on the Lindblad torque corresponds to the angular momentum carried 
H planet's ability to clear a gap around its orbit. This ability varies ^^^y ^y the spiral density waves the planet generates in the disc, 
.^ , with the planet-to-primary mass ratio, the disc aspect ratio (that ^here the flow relative to the planet becomes approximately su- 
is, the disc pressure scale height-to-radius ratio), aiid ttie disc's pgrsonic. In discs with decreasing temperature profiles, it is a 
turbulent viscosity (for more details, see I Crida et al. ' 2006). The negative quantity (iTanakaet alJ2002l:lPaardekooperet al.l2010h . 
migration of planets that do not open a gap is referred to as ^^ich, by itself, would drive type I migration on timescales 
type I migration. It typically applies to low-mass planets, up shorter than a few xlO^ yrs (e.g.,mrd 1997). Local variations 
to a few Earth masses if the central object has a solar mass. In j^ j^e disc's temperatur e and/or density profiles, due for example 
massive protoplanetary discs, intermediate-mass planets build- ^^ o pacity transitions (Menou & Goodman 2004) or to dust heat- 
ing up a partial gap (e.^., Saturn-hke planets) may undergo rapid j^g fflasegawa & Pudritzll2010h mav however change the sign 
runaway migration dMasset & Papaloizou. , 2003 . ) . More massive ^^^j magnitude of the differential Lindblad torque. 
)lanets open a clean gap and are subiect to type II migration 

iLin & P apaIoizoul ll986 l: ICrida & MorbidelUll2007h . The corotation torque accounts for the exchange of angu- 

The reproduction of the mass-period diagram of known lar momentum between the planet and the coorbital flow. It has 

exoplanets by models of planet population synthesis is very two components. One scales with the gradient of disc vortensity 

sensitive to the expression for the tidal torque driving type (vorticity to surface density ratio, projected along the vertical di- 
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rection) across the horseshoe region (I WardI 1 99Tt lMassetll200 ih . 
The other features the entropy gradient across the horseshoe 
regio n ([Baruteau & Masset 2 08; Paardekooper & Pap aloizou 
120081: iMasset & CasolJI 120101: iPaardekooner et alJ 120101 . The 
magnitude of the corotation torque depends on the timescales 
for viscous and thermal diffusion inside the horseshoe region, re- 
quired to maintain the local gradients of vortensity and entropy. 
In the absence of diffusion processes, both gradients would be 
progressively flattened out by advection inside the horseshoe re- 
gion, and the corotation torque would ultimately vanish. This 
is called saturation of the corotation torque. It may be avoided 
by requiring that the viscous and thermal diffusion timescales 
across the horsesh oe region be shorter than the horseshoe li- 
brati on timescale (iPaardekooper et all 120111: iMasset & Casolil 
12010 ). When both diffusion timescales exceed the horseshoe U- 
tum timescale (a fraction of the libration timescale), the coro- 
tation torque is a non-linear process also known as horseshoe 
drag. When they become comparable to, or shorter than the U- 
tum timescale, the corotatio n torque decreases towards its valu e 
predicted by linear theory (IPaardekooper & Papaloizou! l2009h . 
referred to as the linear corotation torque. 

Depending on the disc gradients and diffusion processes at 
the planet orbit, the amplitude of the corotation torque can be 
comparable to, or larger than that of the differential Lindblad 
torque. The impact of the corotation torque on type 1 migration 
can therefore be of considerable importance. A lot of efforts have 
been recently put forward to establish simple and accurate ex- 
pressions for the corotation torque as well as for the differential 
Lin dblad torque, which can be used in population syn thesis mod- 
els (^Masset & Casolj 20 10: Paar dekooper et al.l 2()ll). These ex- 
pressions were obtained for 2D non-magnetic viscous disc mod- 
els. 

The tidal torque has also been investigated in magnetized 
non-turbulent discs. The case of a 2D disc with a toroidal mag- 
netic field has been studied by TerquemI (12003 ) through a lin- 
ear analysis. She found that the differential Lindblad torque is 
reduced with respect to non-magnetized discs (waves propa- 
gate outside the Lindblad resonances at the magneto-sonic speed 
rather than the sound speed), and that there is no linear coro- 
tation torque. Instead, angular momentum is taken away from 
the planet by slow MHD waves propagating in a narrow annu- 
lus near magnetic res onances. These results were confirmed by 
iFromang et al.l (l2005h with non-linear 2D MHD simulations in 
the regime of strong magnetic field (the plasma y6-parameter was 
taken equal to 2 in their study, where yS is the ratio of the thermal 
to magnetic pressures). More recently, 2D and 3D disc models 
with a poloidal magnetic field were investigated bv iMuto et al.l 
( 120081) in the shearing sheet approximation. 

It is unclear how turbulence affects the tidal torque. So 
far, only a couple o f stud ies have investigated this issue. 
iNelson & Papaloizoul (|2004|) performed 3D simulations of lo- 
cally isothermal discs fully invaded by MHD turbulence due 
to the MRI. They found that the running time averaged tidal 
torque on a fixed protoplanet experiences rather large-amplitude 
oscillations, and its final value quite substantially differs from 
the torque value expected in laminar discs. Similar results were 
obtained by Nelson (2005), who allowed the planet orbit to 
evolve. A primary reason for the observed difference between 
the laminar torque and the time-averaged turbulent torque is 
that the 3D MHD sim ulations were not converged in time. 
iBaru teau & LinI (1201 Oh considered 2D isothermal discs subject 
to stochast ic forcing, u sing the turbulence model originally de- 
veloped bylLaHiEniiiral- (2004). They showed that, when time- 
averaged over a sufficient long time period, both the differential 



Lindblad torque and the corotation torque behave very similarly 
as in laminar viscous disc models. 'U ribe et aP (1201 ih have re- 
cently performed 3D MHD simulations of planet migration in 
weakly magnetized turbulent discs, including vertical stratifica- 
tion and a range of planet masses. For type 1-migrating planets, 
they obtained a positive running time-averaged torque in disc 
models where the density near the planet's orbital radius be- 
comes a slightly increasing function of radius due to the disc 
structuring by turbulence, suggesting the existence of an unsatu- 
rated corotation torque in such disc models. 

In this paper, we revisit the properties of the tidal torque 
experienced by a fixed planet in disc models with hydromag- 
netic turbulence generated by the non-linear development of the 
MRI. Locally isothermal discs without vertical stratification are 
considered, and non-ideal MHD effects are neglected. 3D MHD 
simulations using two different codes have been carried out, us- 
ing initial values of the plasma yS-parameter larger than or equal 
to 50. The physical model and numerical setup used in our sim- 
ulations are described in Sect. |2] In Sect. [3] we present results 
of simulations with an inner disc cavity. In viscous disc models, 
an increasing density profile yields a large positive corotation 
torque, which may stall type I migra tion. This is often referred 
to as a planet trap (e.g.. IMasset et al.|[2006b ). Our results of tur- 
bulent simulations show a positive time-averaged tidal torque 
experienced by a planet fixed between the two edges of a cav- 
ity, suggesting the existence of an unsaturated corotation torque 
in MHD-turbulent discs. This existence is confirmed in disc 
models with power-law density profiles, which are described in 
Sect. |4| The time-averaged tidal torque obtained in such discs 
is in decent agreement with the torque value predicted in 2D 
non-magnetized viscous disc models. The discussion in Sect. |5] 
gives more insight into the existence of horseshoe dynamics in 
weakly magnetized turbulent discs, as well as a detailed compar- 
ison between the torque density distributions obtained in MHD- 
turbulent and laminar disc models. This comparison indicates 
the existence of an additional corotation torque in discs with 
a weak toroidal magnetic field, which will be investigated in a 
forthcoming study. Conclusions and future directions are drawn 
in Sect. |6] 



2. Model description 

We explore the properties of the tidal torque between a planet 
and its nascent protoplanetary disc, wherein turbulence is driven 
by the non-linear development of the MRI. A question that needs 
to be addressed is whether an unsaturated corotation torque ex- 
ists and can be maintained in the long term in such turbulent 
discs. For this purpose, MHD simulations have been carried out 
using two different codes, which are described in Sect. 12.1! The 
physical model and numerical setup used in our simulations are 
detailed in Sect. 12.2! The parameters common to all simulations 
presented hereafter are listed in Table [1] 

2.1. Codes 

We investigate the nature of the corotation torque, a process oc- 
curring on the width of the planet's horseshoe region that, for 
planets subject to type I migration, is smaller than the pressure 
scale height (e.g.. IMasset et al.ll2006al) . At this scale, the turbu- 
lence properties can potentially be influenced by details of the 
numerical scheme. Thus, we use two different codes that have 
different dissipation properties. This will help constrain the ro- 
bustness of our results, and their limitations. 
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Table 1. Disc, grid and planet parameters common to Sect. [3] 
(disc with an inner cavity) and Sect. |4] (discs with a power-law 
density profile). 



Parameter 



Value 



Radial resolution 

Azimuthal resolution 

Vertical resolution 

Initial magnetic field 

Aspect ratio 

Planet mass 

Planet location 

Planet's softening length 

Frame 

Radial boundary condition 



320 cells in i? 6 [1;8] 
480 cells in <^ 6 [0;n] 
40 cells in ze [-0.3; 0.3] 
Toroidal, set up in i? 6 [1.5; 5] 
H/R = 0.1 alR = 3 



Mp = 3 X 10-* M,^ 







p - 3, i^p = n/2, Zp 
E = 0.2H(R^) 
Corotating with R = Rp 
Wave-killing zones in 
Re [l;1.5]andi?e [7; 8] 



The two codes that we use to perform the 3D MHD simu- 
lations presented in this paper are NIRV ANA (Ziegler & Yorke 
[1997) and RAMSES JTevssied [20021: Fromang et all 1200^ 
NIRVANA uses an algorithm very si milar to the ZEUS cod e 
to solve the equations of ideal MHD dStone & Norman|[T992[) . 
RAMSES is a finite-volume code that uses the MUSCL-Hancok 
Godunov scheme. For the purpose of this project, a uniform grid 
version of the code (i.e., without any Adaptive Mesh Refinement 
capabilities) was extended to allow the use of cylindrical coor- 
dinates. For comparison, results of 2D hydrodynamical simu- 
lations are also presented, performed with the codes RAMSES 
and FARGO (Masset 2000a b). (None of these 2D simulations 
includes a magnetic field). FARGO is very similar to the 2D 
hydrodynamical version of NIRVANA. Its specificity is to use 
a change of rotating frame on each ring of the grid, which in- 
creases the timestep significantly. 

Results of simulations are expressed in the following units: 
the mass unit is the mass of the central star (M*), the length unit 
is the planet's (fixed) orbital separation (^p), and the time unit is 
the planet's orbital period (Toih) divided by 2;7r. Whenever time 
is expressed in orbits, it refers to the orbital period at the planet 
location. 



2.2. Physical model and numerical setup 

2.2.1. Gas disc model 

We adopt a simple 3D magnetized disc model in which non- 
ideal MHD effects, self-gravity and vertical stratification are ne- 
glected. No explicit kinematic viscosity is included. The ideal 
MHD equations are solved in a cylindrical coordinate system 
{R,(p,z] centred on the central star, with R e [1:8], (p e [0:7r] 
and z e [-0.3:0.3]. The frame rotates with angular frequency 
equal to the Keplerian frequency at R - 3 (the fixed location of 
the planet, when included), and the indirect term that accounts 
for the acceleration of the central star by the planet is included 
in the equations of motion for those simulations that include a 
planet. For simplicity, a locally isothermal equation of state is 
used, with the gas pressure f gas and mass volume density p sat- 
isfying Pgas = pel, where the gas sound speed Cj is specified as a 
fixed function of radius. It is related to the disc's pressure scale 
height H through H = Cs/Qk, with ^k the local Keplerian angu- 
lar velocity. The disc aspect ratio, h = H/R, is defined to be 0. 1 
at R = 3. The vertical extent of the disc model is thus equal to 
two pressure scale heights at the planet's orbital radius. Ideally 
we would prefer to choose a smaller disc thickness, more typical 
of protoplanetary disc models. But the resolution requirement of 



maintaining a turbulent flow driven by the MRI, combined with 
the need to run simulations for many hundreds of planet orbits, 
means that we are forced to adopt a relatively thick disc model. 

The disc is set up in radial equilibrium, with the centrifugal 
acceleration and the radial acceleration related to the pressure 
gradient balancing the gravitational acceleration due to the cen- 
tral star (function of R only). A small level of random noise (5% 
of the local sound speed) is added to each component of the gas 
velocity. The initial radial profiles of the disc's density and tem- 
perature will be specified in Sects.[3]and[4| The initial magnetic 
field is purely toroidal with a net azimuthal flux. Its radial profile 
is set by the choice for the plasma /^-parameter, which is taken to 
be constant throughout the disc region in which the initial field is 
introduced. Several values of/3 in the range [50-400] are consid- 
ered throughout this study, depending on the model. Their value 
will be specified below. The initial magnetic field is introduced 
everywhere except near the radial boundaries, where the field is 
set to zero in /? E [1; 1.5] and R e [5:8]. 

The long term action of the turbulent stresses causes the 
disc density profiles to evolve significantly, leading eventually 
to most of the mass being relocated to the close vicinity of the 
inner and outer radial boundaries. To avoid this and obtain a 
st eady-state density profile, we adopt an approach similar to that 
of iNelson & Gressel (1201 Ol) . and reinforce the original density 
profile on some specified time scale. This procedure, which we 
call the mass-adding procedure for future reference, is achieved 
by solving 



dp 
dt 



(p)^ - Pa 



(1) 



alongside the ideal MHD equations. In Eq. ([T]), (p)^ denotes the 
azimuthally-averaged density (function ofR and z), po is the ini- 
tial density (function of R only), and t is set to 20 local orbital 
periods. We checked by means of 2D hydrodynamical viscous 
simulations with a kinematic viscosity corresponding to a similar 
amount of turbulence as in our MHD simulations that the mass- 
adding procedure has a very small impact on the time evolution 
of the tidal torque, since the planet we consider (see Sect. 12.2.2b 
does not open a gap over the duration of our simulations. Some 
additional remarks on the use of the mass-adding procedure will 
be made in Sect. 



2.2.2. Planet parameters 

We adopt a planet of mass Mp = 3 x lO^^M*, which is held 
on a fixed circular orbit at R^ - 1>, z^ - Q (here and in the fol- 
lowing, all quantities with the subscript p are meant to be evalu- 
ated at the planet's orbital radius). This coiTesponds to a Saturn- 
mass planet if the central object is a Sun-like star. Given the 
large disc aspect ratio in our model (h^ - 0.1), the dimension- 
less parameter determining the flow non-linearity near the planet 
Mp/{hpMi,) = 0.3 < 1 . In addition, the large turbulent stress gen- 
erated by MHD turbulence (the effective viscosity parameter a 
is a few percent, as we shall see in the next sections) prevents 
the planet from opening a gap over the duration of our calcu- 
lations (several hundred orbits). The planet mass we consider 
is therefore relevant to the type I migration regime, although it 
is somewhat larger than would normally be adopted in studies 
of the migration of low-mass planets. The need to obtain con- 
verged estimates of the torque experienced by the planet over 
runs times of several hundred planet orbits, however, has led us 
to adopt this larger than prefeiTed value. The gravitational po- 
tential of the planet is independent of z, this being consistent 
with our cylindrical disc setup. It is smoothed over a softening 
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Fig. 1. Results of 2D viscous disc models with FARGO, using a 
constant kinematic viscosity. The steady running time-averaged 
specific torque on the planet (code units) is plotted against the al- 
pha viscosity at the planet location. Results are obtained for sev- 
eral values of the grid's azimuthal extent Aip but with the same 
resolution. 



length, s, equal to Q.2H{Rp). The calculation of the force exerted 
by the disc on the planet excludes gas located within the planet's 
Hill radius {Ru, the excluded volume is a cylinder of radius Rh 
centred around the planet's position). 



2.2.3. Grid resolution and boundary conditions 

Resolution — The grid resolution in all the simulations pre- 
sented in Sects. [3] and a is (Nr,N^,N,) = (320,480,40). At 
the planet location, the pressure scale height is thus resolved 
by about 15 cells along the radial and azimuthal directions, 
and by 20 cells along the vertical direction. Since we aim to 
investigate the properties of the corotation torque, the radial 
extent of the planet's horseshoe region needs to be sufficiently 
resolved. For type I migrating objects, its half-width, Xs, is 
^ 1.2R^yjMp/(hpM^,) (e.g.. iMasset et al.l l200 6a). Here, it is 
resolved by about 9 cells along the radial direction, which is 
similar to th e resolution used in re cent studies of the corotation 
torque (e.g.. lBaruteau & Linll2010l) . 



Azimuthal domain — As akeady mentioned in Sect. 12. 2. Tl our 3D 
MHD simulations cover an azimuthal extent A(f - n. This choice 
is the best compromise between several limitations. Having 
A(/? = 271 is a natural choice for studying the corotation torque, 
and to our knowledge all works on the corotation torque have 
used this azimuthal extent. But, at fixed resolution, this choice 
requires a maximally large A^^, the number of grid cells along 
the azimuthal direction. At the same resolution, decreasing Atf 
reduces A^^, but at the cost of decreasing the corotation torque. 
This point is illustrated in Fig. [T]for 2D hydrodynamical sim- 
ulations with FARGO. These simulations use the same setup 
as described above for the 3D MHD simulations, except that a 
constant kinematic viscosity v is included. The initial surface 
density and temperature profiles decrease as /? '^^ and 7? ', re- 
spectively. Fig. [U displays the (stationary) running time-average 



of the specific torquqj obtained for three series of runs with 
A<^ = nil, n and In, for which A^^ - 240, 480 and 960, respec- 
tively. For each series, the alpha viscosity parameter at the planet 
location, ap = v/(//p Qk,p), ranges from 6 x 10"^ to 6 x 10"^. At 
a given ap, the corotation torque can be estimated as the differ- 
ence between the total torque, and the total torque at smallest 
dp (which is a good approximation for the differential Lindblad 
torque; note that its value remains almost unchanged by reduc- 
ing Aifi from 2n to njT). The corotation torque is approximately 
maximum when the viscous diffusion timescale across the horse- 
shoe region, TvIsc ~ -^s /^' i^ shorter than the libration timescale. 



Tiib(A9j) = 



8;7r7?p Atfj 

L_ X 1. 

3npjCs 27r ' 



(2) 



but exceeds the U-turn timescale dBaruteau & MasseBl2008l) . 

TU-tuin ~ Zip X Tiib(A<^ = l-K). (3) 

We note that, for A(^ - 2n, the value of a^ at which the corotation 
torque is maximum is in good agreement with the expression 
0.16 (Mp/M*)3''2 /j^4 Qf Baruteau & Lin (2010), which equals ~ 

8 X lO""* for our setup. Decreasing Aip reduces Tub in comparison 
to Tu-turn, which implies that 

- a proportionally larger viscosity is required to maintain the 
same saturation level of the corotation torque (that is, the 
same ratio Tvisc/Tiib), 

- the maximum amplitude of the corotation torque decreases 
as the ratio riib/Tu-turn approaches unity. The linear coro- 
tation torque (that is, the corotation torque in the limit of 
very high viscosities) is not significantly modified by reduc- 
ing A<f. 

Taking Aip = tt is found to be a good compromise between a 
large enough corotation torque, easily measurable in turbulent 
simulations, and a tractable numerical resolution. For compar- 
ison purposes, the 2D hydrodynamical simulations presented 
in Sects. |3] and |4] also have Aip - n. In the following, time 
is expressed in planet orbital periods. Despite the choice for 
A(f = 7T, the planet's orbital period is implicitly calculated as 
27r/Qp. 

Boundary conditions — Periodic conditions are adopted at the 
azimuthal and vertical boundaries. The treatment of the radial 
boundary condition differs in the two codes. Hydrodynamic vari- 
able fluctuations are d amped i dentic ally, using the scheme de- 
scribed in Ide Val-Borro et al.l (|2P06"). This scheme, which re- 
laxes the density and velocity components towards their initial 
value, is applied over R e [1; 1.5] and R e [7; 8]. In order to 
avoid the magnetic field reaching large values near the bound- 
aries, different strategies are used in both codes. In NIRVANA, 
a reflecting boundary condition is applied to the gas radial ve- 
locity that prevents the magnetic field from penetrating into the 
ghost zones. In RAMSES, a magnetic resistivity (?/) is used in 
the radial buffer zones (riin - 5 x 10"^ and ?7out - 10"^ at the 
inner and outer boundaries, respectively). However, such an effi- 
cient field diffusion simultaneously causes the toroidal magnetic 
flux to escape from the grid, ultimately resulting in the volume- 
averaged a value decreasing with time. To avoid this effect, at 
each timestep, an additional axisymmetric toroidal field is added 
in each grid cell to guarantee that the total initial toroidal flux 



' As is conventional, we refer to torque as the torque exerted by the 
disc on the planet. If T denotes the torque, the running time-averaged 
torque, F, is calculated as Fit) = ' ' x f r{u)du. 
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Fig. 2. Surface density as a function of radius shown at different 
times during the simulation with an inner disc cavity. 

is conserved. In doing so, the additional toroidal field is radially 
distributed such that the ratio of its associated magnetic pressure 
and of the initial thermal pressure is uniform over the computa- 
tional domain. 



3. Disc model with an inner cavity 

In locally isothermal discs, the corotation torque is expected 
to be particularly strong in regions where the density gradi- 
ent is large, and for this reason we have performed a sequence 
of simulations with the planet located in a disc with an inner 
cavity. The planet itself orbits in the transition region between 
the outer high-density disc and the inner l ow-density c avity - 
a region often referred to as a planet trap (iMasset et al.ll2006CT 
iMorbidelli et"ani2008l) . 

The initial disc model in this case was set up by first running 
a laminar disc model with initial density and temperature profiles 
scaling as /? '''^ and /? ^ respectively. The functional form for 
the kinematic viscosity, v, used to establish the cavity is 



Vin if R< Rin 

(Vout - Vm)(R - Ria) 



V = Vin + 
V 



{Rom - ^in) 
Vout if ^ > rout, 

1-3 



if Ri„<R< Rn 



(4) 



with Vin = 2 X 10" , Vout - Vin/7, ^in = 2.5 and Rout = 3.5. 
All other parameters (e.g., grid resolution, boundary conditions) 
are as described in Sect. |2] Note that the grid's azimuthal extent 
equals n in all the simulations presented from this section on- 
wards. The numerical simulations presented in the present sec- 
tion were carried out with NIRVANA. 

The laminar simulation was run for -10"^ orbits, after which 
time the surface density profile had achieved the near steady- 
state shown by the smooth (green) curve in Fig. |2] (labelled 
as "Initial profile"). This profile, and the associated azimuthal 
velocity profile, were used to construct a 3D cylindrical disc 
model into which a purely toroidal magnetic field withyS =100 
was introduced throughout the disc. This magnetized disc model 
was evolved for 76 orbits, until turbulence was fully developed 
throughout. The mass-adding procedure described in Sect. 12.2.11 



was used to relax the density profile towards the smooth initial 
density profile throughout the turbulent run. 

The turbulent simulation was restarted by including a planet 
of mass Mp = 3 x 10"'* M* atR - 3. The disc-planet system was 
evolved for ~ 420 orbits. The evolution of the surface density 
profile (vertically-integrated density) during this time is shown 
in Fig. |2] Although there is a modest change in the surface den- 
sity gradient at the cavity edge, due to the MHD turbulence, we 
see that a steep gradient is maintained. The time evolution of 
the volume-averaged contributions to the effective turbulent al- 
pha viscosity (stress to thermal pressure ratios) is shown in the 
left panel of Fig. [3] and it is clear that these globally averaged 
stresses maintain a near steady-state. Snapshots of the variation 
of the total stress as a function of radius are shown in the right 
panel of Fig. |3] and again a near steady-state is observed. The 
effective a profile shows a maximum within the cavity, due to 
the low density and pressure there. A contour plot representing 
the disc midplane density after 127 planet orbits is displayed in 
Fig. |4] showing the position of the planet sitting at the edge of 
the cavity. 

The running time-average of the torque for the turbulent disc 
run is shown by the red curve in Fig. |5] It remains positive, 
clearly indicating that the corotation torque is strong and main- 
tained in an unsaturated state throughout the simulation. For 
comparison, we ran four additional 2D simulations using lam- 
inar disc models with different values of the kinematic viscosity. 
The initial surface density profile is the same as the initial profile 
shown by the green curve in Fig.|2] For comparison, the torques 
obtained in the 2D simulations have been multiplied by the scal- 
ing factor / = 0.6 to account for the vertical extent of the 3D 
simulations, where z extends from -0.3 to 0.3. We see that the 
laminar disc model with a^ - 0.01 (shown by the green curve in 
Fig. |5]l maintains a total torque whose cumulative time-average 
is similar to that of the turbulent run, but note that the surface 
density profiles in both runs differ about the planet location, as 
can be seen in Fig.|2] 

The laminar run represented by the blue line in Fig. |5] 
adopted the kinematic viscosity profile given by Eq. (|5]l, such 
that V = 1.14 X 10"^ at the planet location, corresponding to 
the rather large value of a^ ~ 0.6. The maximum unsaturated 
corotation torque is expected to occur when o-p ~ 2 x 10"^ (see 
Fig.[TJ, and for values larger than this, the corotation torque ap- 
proaches a smaller value predicted by linear theory. Our results 
are consistent with this. The magenta line in Fig. |5] shows the 
result from a laminar simulation with o-p - lO"-', a value for 
which we expect the corotation torque to be partially saturated. 
This expectation in confirmed by the simulations, and we see 
that the torque obtained with ffp = 10"^ lies between those for 
Qfp = and Op = 10"^. The black curve plotted in Fig.|5]is the 
torque obtained for an inviscid laminar disc. Clearly the net neg- 
ative torque demonstrates that the corotation torque in this case 
becomes fully saturated, leaving only the differential Lindblad 
torque. 

These cavity runs show that MHD turbulence with o-p ^ 0.02 
allows the corotation torque to be sustained at a value close to its 
maximum, unsaturated value for the disc and planet parameters 
that we have considered. Significant changes in the effective a, 
due to the turbulence being more or less vigorous, however, will 
cause the corotation torque to weaken. Our runs demonstrate that 
a planet trap maintained in a turbulent protoplanetary disc can 
be effective in preventing the large scale migration of embedded 
protoplanets. 
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Fig. 3. Left: volume-averaged alpha values for the disc model with an inner cavity, shown from the point in time when the planet 
is inserted in the disc. As in all other panels, time is expressed in orbital periods at R - 3 (the planet's fixed location). Right; 
azknuthally- and vertically-averaged profiles of alpha for the same disc model, shown at different times after the planet inclusion. 
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As a conclusion to this section, our results strongly suggest 
the existence of an unsaturated corotation torque in turbulent 
protoplanetary discs with weak magnetic field (here^S = 100). 
Although the agreement with 2D viscous laminar simulations is 
quite good, the complexity of the disc structure resulting from 
the presence and time-variation of the cavity makes it difficult to 
perform a quantitative analysis. The strongly increasing density 
profile about the planet location, the precise slope of which sig- 
nificantly impacts the corotation torque, and the fact the density 
profile can be approximated as a power-law function of R only 
over a limited region (< +3H) about the planet location hinder a 
quantitative comparison between the torques displayed in Fig.|5] 



and analytic torque formulae. In Sect.|4l a suite of idealized sim- 
ulations using power-law density profiles are constructed that 
will help study quantitatively the corotation torque properties. 

4. Disc models with a power-law density profile 

The results of the turbulent disc model with an inner cavity pre- 
sented in Sect. |3] strongly suggest the existence of an unsatu- 
rated corotation torque in discs with hydromagnetic turbulence 
generated by the MRI. In the present section, we investigate 
the corotation torque properties in turbulent discs with initial 
density profiles that are power-law functions of radius. In con- 
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Fig. 5. Running time-averaged torque obtained with an inner 
disc cavity: comparison of laminar and turbulent disc models. 

Table 2. Disc parameters of the turbulent runs in Sect. |4] where 
initial power-law density profiles are adopted. Results obtained 
with Model 1 and Model 2 are presented in Sect. 14.11 and 14.21 
respectively. 



Parameter 



Value 



Initial gas density po '^ R '''^ (Model 1) 

po 0^ R'^'^ (Model 2) 
Initial gas temperature Tq oc /?"' (Model 1) 

To uniform (Model 2) 
Plasma /3-parameter j3 = 50 (NIRVANA), /? = 400 (RAMSES) 



trast to Sect. [3] the simulations in this section have been car- 
ried out with both NIRVANA and RAMSES. A constant initial 
plasma jS-parameter is adopted throughout the disc, which is 
equal to 50 in the runs performed with NIRVANA, and to 400 
with RAMSES. As will be shown below, these different values 
are found to give similar volume-averaged alpha values in both 
codes. 

Two disc models are considered in this section, whose pa- 
rameters are summarized in Table |2] In Model 1 , the initial den- 
sity and temperature profiles, denoted by po and To, decrease as 
R^^^^ and /?"', respectively. In Model 2, po ^ R^^^^ and Tq is 
uniform. In both models, we use pq(R = 1) = 10""* in code units. 
Assuming a Sun-like star, and that the star-planet separation is 
5 AU, this corresponds to an initial surface density at the planet 
location ^ 180 g cm"^ in Model 1, and ^ 60 g cm"^ in Model 2. 
All other parameters are those described in Sect. [2] Recall in par- 
ticular that Tq(R - Rp) is chosen such that hj, = 0.1. The results 
of both models are described in Sects. I4n and l4. 21 respectively. 
In the absence of magnetic field and turbulence, the corotation 
torque would vanish in Model 2, but not in Model 1 . However, 
we do not know whether and how magnetic effects modify the 
dependence of the corotation torque with density and tempera- 
ture gradients. The tidal torque values obtained in both models 
can at best hint whether MHD turbulence impacts the differen- 
tial Lindblad torque and the corotation torque, but they cannot 
allow us to work out quantitatively how each torque is altered by 
MHD turbulence. 



4.1. Model 1: initial density oz R ^^^ and temperature oc R ■ 

The initial density and temperature profiles in Model 1 are the 
same as those of the 2D simulations of Sect. 12.2.31 where the 
influence of the grid's azimuthal extent on the amplitude of 
the corotation torque is examined. The time evolution of the 
volume-averaged alpha parameter, associated with the sum of 
the Reynolds and Maxwell stresses, is shown in the left panel 
of Fig. |6] It is denoted by (fftot)- The point in time from which 
the planet is inserted is shown by a vertical dashed line for each 
simulation. The time to reach a turbulent saturated state is longer 
in the RAMSES run, since it uses a larger/? value. In the satu- 
rated state, (atot) ~ 0.02 in both the NIRVANA and RAMSES 
calculations. At first, this agreement might seem surprising since 
the initial toroidal flux is diff'erent in the two simulations. This 
flux is indeed co nserved throu ghout the simulations, and it is 
well known (e.g.. ' Nelsonll2005h fliat, in the saturated state of the 
turbulence, (o-tot) scales with the initial toroidal flux. This agree- 
ment is in fact due to the different algorithm used by the two 
codes. Being less diff'usive, Goduno v codes are known to give 
in general larger alpha values (e.g.. ISimon et al.ll2009h . In this 
study, we tuned the respective initial toroidal fields in order to 
obtain the same saturated transport properties. 

The time-averaged radial profiles of the alpha parameters as- 
sociated with the Reynolds, Maxwell and total stresses (averaged 
over (fi and z) are displayed in the right panel of Fig. |6] They are 
respectively denoted by ^Rey, fFiviax and fftot- The time average is 
done between 200 and 300 orbits for the RAMSES run, and be- 
tween 250 and 350 orbits for the NIRVANA run. Note the slight 
decrease in both aRey and ffMax near the planet's orbital radius, 
which accounts for the slight increase in the disc density (or, 
equivalently, the disc thermal pressure) at the same location, as 
shown in Fig. |7] We checked that the time-averaged magnetic 
stress is slightly increased inside the planet's Hill radius and at 
the location of the planet's wake, where the magnetic field tends 
to be compressed and ordered. T his is reminiscent of the findings 
by iNelson & Papaloizoul (l2003h (see their figure 13), who how- 
ever considered a planet ~ 15 times more massive than ours. 
At the planet location, cFtot ~ 0.03. This translates into an aver- 
aged timescale for viscous diffusion across the horseshoe region 
^ Xj (fftotCs//)p' ~ 2.5rorb, which is similar to the horseshoe 
U-turn timescale, Tu-turn - '^■Torb, see Eq. (O. We thus expect 
the corotation torque to adopt a value close to that predicted by 
linear theory. 

The time- and azimuthally-averaged disc midplane density is 
depicted in Fig.|7] The time-average is done over the same time 
interval as for the total alpha parameter shown in the right panel 
of Fig. |6] The mass-adding procedure described in Sect. 12.2.1! 
which relaxes the density profile towards its initial value over 
20 local orbital periods, implies that we get a stationary density 
profile slightly reduced compared to the initial one, depicted by 
a dashed curve in Fig.|7] Note that the decrease in the disc mass 
is larger in the RAMSES simulation at all radii. At the planet 
location, the averaged-to-initial density ratio is x; 0.95 and x; 
0.88 for the NIRVANA and RAMSES calculation, respectively. 

Density contours in the disc midplane obtained with the 
NIRVANA run after 210 orbits are shown in Fig.|8] The planet 
is located at x - Q,y - 3, and the color scale has been adjusted 
such that the maximum density contour is the maximum value of 
the initial gas density. The spiral density waves generated by the 
planet are hardly visible, their density contrast being very similar 
to that of the turbulent fluctuations. Time-averaged density con- 
tours are displayed in the left panel of Fig.|9] The time averaging 
is done over 35 orbits, starting 345 orbits after the restart time 
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Fig. 6. Left: time evolution of the volume-averaged total alpha parameter, (o'tot), obtained with Model 1 . Vertical dashed lines indi- 
cate the point in time when the planet is inserted in each simulation. Right: azimuthally- and vertically-averaged alpha parameters 
associated with the Reynolds stress (dotted curves), the Maxwell stress (dashed curves) and the total stress (solid curves) with the 
same model, time-averaged from 200 and 300 orbits for the simulation with RAMSES, and from 250 to 350 orbits for the simulation 
with NIRVANA. 
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Fig. 7. Azimuthally- and time-averaged midplane density pro- 
file obtained in Model 1 with NIRVANA (black curve) and 
RAMSES (red curve). The time averaging is done between 250 
and 350 orbits for the NIRVANA run, and between 200 and 300 
orbits for the RAMSES run. The dashed curve displays the initial 
density profile common to both simulations. 



(or, 520 orbits after the beginning of the simulation). This time 
interval corresponds to ~ 3 horseshoe libration periods. While 
the grid's full azimuthal range is shown by the y-axis, the x-axis 
displays a narrow range of radii about the planet's orbital radius. 
The planet is located at Rp -3 and ipp = n/l, and its wake is 
clearly visible. Overplotted by solid lines are streamlines in the 
disc midplane, evaluated from the time-averaged values of the 
radial and azimuthal components of the gas velocity field over 
the same time interval. Horseshoe streamlines are clearly seen, 
which delimit a well-defined mean horseshoe region, encom- 
passed by mean circulating streamlines. The use of a quite long 
(35 orbits) time interval for the time averaging is aimed at dis- 
playing smooth averaged streamlines, which helps estimate the 
half- width of the planet's mean horseshoe region, Xs- As will be 
shown in Sect. 15.11 horseshoe streamlines are readily observed 



when averaging over 5 orbits (see Fig.[T2b. The blue streamlines 
depict the approximate location of the separatrices of the mean 
horseshoe region. Its half-width at (fi - cpp = I rad is Xs ~ 0.20. 
For comparison, gas surface density contours obtained with the 
FARGO viscous run are shown in the right panel of Fig. |9] along 
with instantaneous streamlines. We see that the width and shape 
of the mean horseshoe region in the MHD turbulent simulation 
are very similar to their counterpart in the laminar viscous simu- 
lation. In the viscous run, the half- width of the horseshoe region 
is Xs ~ 0.23 at(p-(pp - 1 rad. This is in good agreement with the 
value obtained from the expression given bv lPaardekooper et al.l 
(2010) (their Eq. 44), 



Xs = l.lRp X 



0.4 



1/4 



1/2 



fl -0.215. 



(5) 



Fig. [To] depicts the running time-averaged torque (solid 
curves) obtained with NIRVANA (left panel) and RAMSES 
(right panel). We briefly recall that all torques shown in this pa- 
per refer to the specific torque exerted by the disc on the planet, 
and that the content of the planet's Hill radius is excluded in 
the torque calculation. In both runs, the torque reaches a sta- 
tionary value after typically 100 to 200 orbits. Although not 
shown here, this timescale is consistent with the amplitude of 
the stochastic torque standard deviation being similar to that of 
the ty pe I migration torque in an e quivalent lamina r disc model 
(e.g.. (Nelson & Papaloizoul 120041: iBaruteau & Linll2010l) . This 
similarity can be qualitatively expected from the similar density 
contrasts of the turbulent fluctuations and of the planet wake, 
illustrated in Fig. [8] We comment that the stationary running 
time-averaged torque obtained with both codes are in very good 
agreement, that of the RAMSES run being 10 to 20% smaller 
than that of the NIRVANA run. This slight difference may be (at 
least partly) attributed to the time-averaged density profile being 
slightly smaller in the RAMSES run (see Fig.|7]i. 

How do the averaged torques of our MHD turbulent simula- 
tions compare with the predictions of laminar disc models? To 
address this question, we have performed two series of 2D lami- 
nar disc calculations using FARGO and RAMSES. In one series, 
the disc is inviscid. In the other series, a viscous disc model is 
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Fig. 8. Density contours in the disc midplane shown at 210 orbits for the NIRVANA run with Model 1. The planet is located at 

x = Q,y^3. 



c 



: 



4-1 0'= 5-10"'' 6-1 0"' 7-1 0"^ 8-1 0"'^ 9-1 0"' MO" 4-1 0"^ 5-1 0"= 6-1 0"° 7-10^ 8-1 0"'* 9-1 0"^ MO" 

Density averaged from 520 to 555 orbits Surface density at t = 520Torb 





Fig. 9. Left: contours of the midplane density obtained in Model 1 with the NIRVANA run, and time-averaged between 520 and 
555 orbits. Streamlines in the disc midplane, time-averaged over the same time interval, are overplotted by solid curves. Right: gas 
surface density obtained at 520 orbits with disc Model 1 in the FARGO viscous run, with instantaneous streamlines overplotted. In 
both panels, the horizontal dashed line shows the planet's corotation radius, defined as the location R^ where Q(^c) - ^p- 



used with a prescribed kinematic viscosity v. In the latter case, 
the FARGO run uses a fixed radial profile for a that matches 
the time-averaged total alpha profile obtained in the MHD run 
with NIRVANA (black curve in the right panel of Fig. |6ll. The 
RAMSES run uses a constant kinematic viscosity, chosen such 
that ffp equals the averaged alpha value obtained in the RAMSES 



MHD run at the planet location (a 0.03). In both series of lam- 
inar disc models, the initial density profile is the time-averaged 
density profile of the corresponding MHD simulations (shown 
in Fig.|7]). We have found that the steady-state density profiles in 
the viscous simulations are reduced by about 10% with respect 
to their initial profile. At the end of the inviscid simulations (400 



C. Baruteau et al.: Corotation torque in weakly magnetized turbulent discs 



3 

E 

o 
4-« 


-2.0.1 0-' 
-4.0.10' 


■a 
m 


-6.0.10' 


9 


-8.0.10' 


E 


-1.0.10"' 


•*-» 




c 


-1.2.10"' 


c 




1. 

3 


-1.4.10"' 




-1.6.10' 





g, -6.0-10"' 



> 
ni 

I -1.0-10"' 

g' -1.2-10' 

"E 

i -1.4-10"' 

DC 



-1.6-10 



Laminar disc with a = 0.03 



Laminar disc witli a = 



100 200 300 400 500 
Time (orbits after restart) 



50 100 150 200 250 300 
Time (orbits after restart) 



Fig. 10. Running time-averaged specific torque obtained with Model 1 , for the NIRVANA run (left panel), and the RAMSES run 
(right panel). For comparison, the final value of the time-averaged torques obtained from 2D laminar disc models are overplotted, 
wherein the initial density profile is the time-averaged density profile obtained with the MHD run. Note that in the left panel, the 
turbulent run with NIRVANA is compared to laminar runs caiTied out with FARGO. Two series of laminar disc models have been 
performed, one with an inviscid disc (ap - 0), and one with a viscous disc (corresponding to Op = 0.03). 



orbits), a tiny gap has opened about the planet's orbital radius, 
with a density contrast ^ 20%. 

The final value of the running time-averaged torque obtained 
in both series of laminar calculations is overplotted in Fig. [10] 
by dashed lines. The 2D laminar torques have been multiplied 
by the scaling factor / - 0.6 to account for the vertical ex- 
tent of the 3D simulations (where z ranges from -0.3 to 0.3). 
The torques obtained in the turbulent and viscous disc models 
are in good agreement, their relative difference being less than 
10% for the NIRVANA and FARGO runs, and less than 20% for 
the RAMSES runs. Taking into account that the stationary den- 
sity profile of the viscous runs is ^ 10% smaller than the time- 
averaged density profile of the turbulent runs, the actual rela- 
tive difference between the viscous and turbulent torques ranges 
from 20% to 30%. In the series of inviscid runs, the torque is 
more negative, as expected, since the (positive) corotation torque 
saturates after a few libration timescales, leaving only the differ- 
ential Lindblad torque. 

It is instructive to compare our results with the predictions o f 
the analytical torque expressions of Paardekoop er et al.l (1201 Oh. 
derived for 2D viscous disc models with power-law density pro- 
files. At the end of the inviscid simulations, although a small 
gap has opened, the density profile can still be approximated as 
a power-law profile. For instance, in the FARGO inviscid run, 
it is fairly well approximated by HpiR/Rpy^^^ with Sp ^ 5.2 x 
10~^. The torque expression in Eq. (14) of iPaardekooper et aP 
(1201 0') predicts a specific differential Lindblad torque Fl - 
-6.8(Mp/M*)/!p^Sp/?pQ^. Multiplying this torque expression 

by the aforementioned scaling factor / gives Fl ~ -1.9 x 10"^ 
for the density profile of the FARGO inviscid run. This pre- 
dicted value compares decently with our simulation result (x 
-1.45 X 10"^, see left panel of Fig.fTOb. Given the large viscosity 
parameter in our turbulent model, we expect that the corotation 
torqu e is close to its linear value (IPaardekooper & Papaloizoul 
I2OO 9I). Similarly, using the parameters of the FARGO viscous 
run, Paardekooper et al. (2010) predict a linear corotation torque 
«; 5.6 X 10"^ (see their Eq. 16). This is in good agreement with 
our findings, where the corotation torque (estimated as the torque 



of the viscous model substracted from the torque of the inviscid 
model) amounts to a 6.5 x 10"^. 

The results of this section brings more support to the exis- 
tence of an unsaturated corotation torque in weakly magnetized 
turbulent discs. We find a decent agreement between the run- 
ning time-averaged tidal torques in this particular MHD turbu- 
lent disc model, and in a viscous non-magnetic disc model with 
similar viscosity parameter at the planet location. The proper- 
ties (shape, width) of the mean horseshoe region compare de- 
cently with those of the non-magnetic case. We cannot conclude, 
however, that the differential Lindblad torque and the corotation 
torque take similar values in both types of models. It is possi- 
ble that MHD turbulence alters both the differential Lindblad 
torque and the corotation torque, in such a way that the total tidal 
torque is hardly changed. The results of Model 2, where laminar 
disc models predict no corotation torque, will help us evaluate 
to what extent MHD turbulence alters the differential Lindblad 
torque, and also whether or not the inclusion of magnetic effects 
can introduce an additional torque. 

4.2. Model 2: initial density oz R^^l^ and uniform temperature 



The results of Model 1, presented in Sect. 14.11 strongly indi- 
cate the existence of an unsaturated corotation torque in weakly 
magnetized turbulent discs. The tidal torque compares decently 
with the value expected in laminar viscous disc models. In this 
section, we consider a disc model with a constant temperature, 
and an initial density profile decreasing as R^^^^. We refer to this 
model as Model 2. Laminar non-magnetic disc models predict 
no corotation torque with such disc profiles, and that the tidal 
torque therefore coincides with the differential Lindblad torque. 
For the alpha values of our turbulent disc models (a few percent), 
the differential Lindblad torque is expected t o have no signifi- 
cant d ependence with the disc viscosity (e.g., iMuto & Inutsukal 
I2OO9I) . 

As in Model 1, we have caiTied out 2D inviscid and viscous 
simulations using the FARGO and RAMSES codes. Their initial 
density profile is the time-averaged profile obtained in the cor- 
responding 3D MHD runs (that of the NIRVANA run is used as 
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Fig. 11. Running time-averaged specific torque obtained in Model 2 with the NIRVANA run (left panel) and the RAMSES run (right 
panel). As in Fig. [TO] the results of 2D inviscid and viscous disc models are overplotted, wherein the initial density profile is the 
time-averaged density profile of the coiTesponding MHD run. 



initial condition for the 2D run with FARGO). The FARGO vis- 
cous run imposes a radial profile of alpha viscosity equal to the 
radial profile cPtot of the NIRVANA run. In the RAMSES viscous 
run, a constant kinematic viscosity is used, such that the corre- 
sponding value of Q-p - 0.03, the value of (Op) in the RAMSES 
MHD simulation. 

The time evolution of the running time-averaged torque is 
depicted in Fig.[n] Results with NIRVANA and FARGO are dis- 
played in the left panel, those with RAMSES in the right panel. 
As in Fig. (TO] the torques of the 2D runs are multiplied by the 
scaling factor / = 0.6, to compare them with the torques of the 
3D runs. The stationary running time-averaged torques of both 
MHD simulations are in close agreement. This is maybe partly 
fortuitous, given that the corresponding time-averaged density 
profiles differ by about 25% (not shown here, but analogous 
to the different behaviors obtained in Model 1, and shown in 
Fig-13- This ^ 25% relative difference in density, however, is 
consistent with the relative difference between the torques in the 
laminar viscous simulations. The tiny torque difference between 
the inviscid and viscous runs with FARGO is most probably due 
to the opening of a shallow gap structure about the planet loca- 
tion. Gap formation also occurs in the RAMSES run, but in this 
particular case, it is found to have a negligible influence on the 
stationary running time-averaged torque. 

For the paramete rs of Model 2, the differential Lindblad 
torque expression of iPaardekooper et al.l (2010) (see their Eq. 
14) reads Fl - -3.9{Mp/M^)hfl.pR*^Q.l. The time-averaged 
density profiles of our turbulent simulations can be approx- 
imated as "LpiR/Rpy^^^ near the planet location, with Ep ^ 
1.9 X 10-5 and 1.6 x 10"^ for the NIRVANA and RAMSES 
runs, respectively. The value of Fl predicted by laminar disc 
models then amounts to about -3.9 x 10"^ and -3.3 x lO"*" for 
the NIRVANA and RAMSES runs, respectively. These predicted 
values agree with our numerical results to within 10% for the 
NIRVANA run, and 30% for the RAMSES run. 

The running time-averaged torque of the NIRVANA tur- 
bulent run is in reasonable agreement with the torque of the 
FARGO viscous run, their relative difference being about 30%. 
The agreement is less good for the RAMSES run, where the 
relative difference is ~ 60%. Interestingly, the averaged torque 
measured in both MHD simulations is more negative than in the 



laminar viscous runs, while the opposite trend is observed in 
Model 1. This comparison suggests that either the differential 
Lindblad torque is slightly increased, or that there exists an ad- 
ditional torque in magnetized turbulent discs. More insight into 
this hypothesis is given in Sect. 15.21 

5. Discussion 

5.1. Horseshoe dynamics 

Although the results presented in previous sections highlight 
the action of unsaturated corotation torques, we have not so far 
demonstrated the existence of horseshoe dynamics in our sim- 
ulations without time-averaging over long timescales. A key is- 
sue is whether the horseshoe trajectories exist and can be de- 
tected easily in the presence of MHD turbulence and the associ- 
ated turbulent velocity fluctuations. A related issue is the ques- 
tion of how rapidly the turbulence causes diffusion across the 
horseshoe region. This is a key factor in determining the mag- 
nitude of the corotation torque. As discussed in Sect. 12.2.31 the 
corotation torque in laminar viscous disc models should have its 
maximal non-linear value when the diffusion timescale across 
the horseshoe region, Tvisc, is shorter than the horseshoe libra- 
tion time, r^b, and is larger than the horseshoe U-turn time 
Tu-turn - hpTin,{A(f = 27t). If the diffusion time Tvisc < Tu-mm, 
however, the corotation torque will approach the lower value ex- 
pected from linear theory. 

In order to examine the diffusion rate across the horseshoe 
region, we have performed simulations in which a small patch 
of the disc in the horseshoe region is polluted initially with a 
passive contaminant that is advected with the flow. We define a 
function tfr(R, if) = C{R, ip) p, where C{R, if) is the concentration 
of the contaminant, and evolve the equation 



^ + V ■ (^v) = 0. 



(6) 



We have carried out three simulations to examine the effects of 
turbulent diffusion. The disc model in each run is the same as 
Model 1 (see Sect. 14. lb . In the first run, no planet was included in 
the simulation. In the second run, a planet with Mp = 3 x 10 "^M* 
was included, and the system was relaxed for approximately 40 
planet orbits prior to the contaminant being introduced. In the 
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Fig. 12. Evolution of the concentration of the passive contaminant, C(R, (p), defined in the text, with superimposed averaged stream- 
lines in the midplane for disc models with parameters equal to those of Model 1 . The contour levels are defined for each panel 
separately to highlight the region occupied by the contaminant at each time. The upper three panels are for a disc with no embedded 
planet. The middle three panels have an embedded planet with Mp = 3 x 10 '*M*. The lower panels are for an embedded planet 
with Mp = 10^^ Mi,. Times are given in units of the planet orbital period. Results are obtained with NIRVANA. 



third simulation, a Jovian mass planet with Mp - 10""' M* was 
included, and a similar process of disc relaxation was undertaken 
to ensure that horseshoe trajectories were fully established prior 
to the release of the contaminant. 

The results of these calculations are shown in Fig. [12] which 
displays contours of the contaminant concentration in the disc 
midplane, at three different times. The top panels show the run 
with no planet included. The initial patch of contaminant is lo- 
cated at R - 3.08, between 1.9 < if < 2.0, and at all heights in 



the disc. Moving from left to right through the panels we see that 
the turbulence causes quite rapid diffusion in the radial direction, 
with the contaminant diffusing over a radial distance ^ 0.2 dur- 
ing a time of 1.41 orbits. The Keplerian shear combines with the 
radial diffusion to cause the contaminant to be stretched along 
the azimuthal direction through advection. Nonetheless, we are 
able to observe that on average the initial patch of contaminant 
continues along its original trajectory, shearing past the nomi- 
nal position of the planet that is indicated by the filled white 
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circle. We have superimposed streamlines calculated from the 
time-average of the velocity field in the disc midplane over 5 or- 
bital periods measured at /? = 3 (the orbital radius of the planet, 
if present). As we are working in the rotating frame centred on 
this radius, these show the Keplerian shear for disc regions that 
are interior and exterior to R = 3, but it is also clear that irregular 
motions associated with the turbulence have not been completely 
averaged out over this time. 

The middle panels of Fig. [12] show the run with Mp - 
3 X lO^^M*. It is immediately obvious that the time-averaged 
velocity field displays the existence of horseshoe streamlines. 
In general, we find that time averaging is required for temporal 
windows of between 3-5 orbits before the horseshoe stream- 
lines become clearly defined in Model 1, for this planet mass. 
The half-width of the horseshoe region is approximately equal 
to 0.2, in agreement with the previous determination from the 
left panel of Fig.|9] where streamlines are time-averaged over 35 
planet orbits. 

We see that the evolution of the passive scalar is rather dif- 
ferent from that observed in the no-planet case, and although 
there is significant radial diffusion we also see clear evidence 
of the fluid following the horseshoe streamlines. In particular, 
the presence of significant contaminant near the downstream in- 
ward separatrix (that is, the region bounded by 2.8 < R < 3 
and 7r/2 < i^ < 2) in the second and third panels demon- 
strates that the fluid has undergone a U-turn after approaching 
the planet. Similarly, at f ~ LSlTorb, the concentration in the 
region bounded by 3.0 < 7? < 3.4 and 0.5 < tp < n/2 is smaller 
than in the no planet case. For our disc model and planet param- 
eters, Tiib - lOroib, Tu-tura - 2rorb, and Fig. [12] indicates that 
Tvisc ^ T"u-tum- Notc that turbulcnt diffusion of a passive scalar 
is powered by the Reynolds stress only. However, the turbulent 
diffusion of vortensity, which should control the magnitude of 
the corotation torque by analogy with viscous disc models, de- 
pends a priori both on the Reynolds and the Maxwell stresses. 
The timescale for turbulent diffusion of vortensity should thus be 
shorter than for the contaminant, and we expect the corotation 
torque to have a value smaller than its fully unsaturated non- 
linear value, in agreement with our earlier discussion of Model 
1. 

The lower panels of Fig. [12] show results for the run with 



M„ 



10 Mi,, and we include this simulation in our discus- 



sion only for illustration purposes. The horseshoe streamlines 
are now much more clearly defined, and the trajectory of the 
passive contaminant follows these streamlines even more clearly 
than in the Mp = 3 x IQ^'^M-t, run. Here, advection of the pas- 
sive contaminant around the horseshoe streamlines occurs more 
rapidly than radial diffusion. For the time scales covered by the 
panels in Fig. [12] the contaminant therefore remains more or less 
contained within the horseshoe region. 

For the disc and planet parameters of our model, we find 
clear evidence of the existence of horseshoe dynamics in weakly 
magnetized turbulent discs, in the sense that mean horseshoe 
streamlines are exhibited along which fluid elements perform 
horseshoe U-turns despite substantial turbulent diffusion. We 
comment that we have not shown however whether it is possible 
or not to get a fully unsaturated horseshoe drag in such turbulent 
configurations. Also, the presence of horseshoe dynamics in tur- 
bulent discs can be attributed to the averaged turbulent velocity 
fluctuations being of smaller amplitude than the U-turn drift rate. 
It is unclear whether or not a horseshoe dynamics still exists for 
those smaller planets for which the U-turn drift rate is smaller 
than the typical amplitude of the gas turbulent velocity fluctua- 



tions. This regime deserves a dedicated investigation that goes 
beyond the scope of this paper 

5.2. Averaged torque density distribution 

The results of simulations presented in Sect.|3]and|4]support the 
existence of an unsaturated corotation torque in weakly mag- 
netized turbulent discs. The running time averaged torques ob- 
tained with the NIRVANA and RAMSES codes are in very good 
agreement overall. Given the small scales we are interested in, 
this is reassuring as to the robustness of our results. 

For the two disc models in Sect. [4] we find a decent agree- 
ment overall between the steady running time-averages of the 
tidal torque in the turbulent and viscous simulations. In Model 
1, the torque relative difference between turbulent and viscous 
models is less than 30% for both codes. In Model 2, however, it 
is less than 60%. To provide more insight into these torque dif- 
ferences, we compare in Fig. [T3] the time-averaged torque den- 
sity profiles obtained for both models (depicted by plus sym- 
bols), with the torque density profiles of the viscous and inviscid 
laminar disc models (star and diamond symbols, respectively). 
Results with FARGO and NIRVANA are shown in the left panel, 
those with RAMSES in the right panel. The torque profiles of the 
MHD runs are time-averaged over 100 orbits. In each plot, the 
approximate location of the separatrices of the planet's horse- 
shoe region is indicated by dashed lines. 

We first compare the torque density distributions obtained in 
the laminar disc models. In Model 1, the torque distributions of 
the viscous and inviscid runs mostly differ within the planet's 
horseshoe region. This difference highlights the positive unsatu- 
rated corotation torque of the viscous disc model. Slight differ- 
ences can be noticed between the corotation torque distributions 
of the FARGO and RAMSES simulations. In Model 2, where 
there is no corotation torque, the torque distributions of the vis- 
cous and inviscid simulations are in very good agreement. Their 
tiny differences show the influence of viscous diffusion on the 
radial profile of the differential Lindblad torque. 

We now compare the averaged torque distribution of the 
MHD runs (plus symbols) with the torque distribution of the vis- 
cous runs (star symbols). Outside the planet's (mean) horseshoe 
region, both distributions are found to be in overall good agree- 
ment. This is a clear demonstration that the differential Lindblad 
torque remains essentially unchanged by the full development 
of MHD turbulence. Inside the planet's horseshoe region, how- 
ever, differences are significant. In Model 2, the turbulent torque 
distribution takes more negative values over most of the horse- 
shoe region, except near the outer separatrix, where it takes more 
positive values. Note the excellent agreement between the turbu- 
lent torque distributions of the NIRVANA and RAMSES runs. 
In contrast, the turbulent torque distribution obtained in Model 
1 tends to take more positive values over the horseshoe region, 
except near the inner separatrix, although this trend is less evi- 
dent in the RAMSES run. These torque differences suggest the 
existence of an additional torque in weakly magnetized turbu- 
lent discs that originates from the planet's horseshoe region. The 
distribution of this additional torque, in particular its asymmetry, 
accounts for the torque differences between the MHD turbulent 
and viscous simulations underlined in Sect. l4.1l and l4.2l 

It is quite surprising that the density distribution of the addi- 
tional torque is not symmetric with respect to the location of 
the planet, or of the planet's corotation radius. This offset is 
obtained with our two codes, but it differs in Models 1 and 2. 
Notwithstanding the offset, both models point out that the ad- 
ditional torque tends to take negative values in the inner half 
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Fig. 13. Time average of the torque density distribution, obtained with Model 1 (top row, plus signs) and Model 2 (bottom row, 
plus signs). For comparison, the torque density distribution of the inviscid and viscous non-magnetic disc models are overplotted 
(diamond and star symbols, respectively). The torque profiles of the turbulent runs have been time-averaged over 100 orbits. In 
each panel, torque density profiles are divided by the maximum of the (absolute value of the) inviscid, viscous, and turbulent torque 
profiles. In x-axis, orbital separation is depicted from R = 2 to R - 4, which corresponds to a radial extent ^ +3//p about the planet 
location. In each panel, the dashed lines show the approximate location of the separatrices of the planet's horseshoe region, and the 
dotted line that of the planet's corotation radius. 



of the horseshoe region, and positive values in the outer half. 
This suggests that the additional torque i s probably not p ow- 
ered by magnetic resonances. As shown bv lTerquemI ( |2003|) . the 
torques at the magnetic and Lindblad resonances have the same 
sign. The torque exerted by the disc on the planet near the inner 
(out er) magnetic resonanc e is therefore positive (negative) (see 
also iFromang et al]|2005l their Fig 6). The torque distribution 
of this so-called magnetic torque, and of the additional torque 
in our magnetized turbulent disc models, therefore have oppo- 
site signs. We also point out that magnetic resonances are barely 
reso lved in our setup . Their location, /?„, as given by linear the- 



ory (lTerquemll20dl . satisfies \R„,-Rj,\ = (2i/p/3)(l + P)-^'^, 
with /3 the effective plasma parameter evaluated from the aver- 
age of the magnetic field squared. For the NIRVANA runs, we 
find j6 X 40, and {R^ - Rp\ ^ 0.03 is scarcely larger than our 
mesh size (x; 0.02). 

It is possible that the additional torque found in our mag- 
netized turbulent disc models is an additional corotation torque. 
Neglecting turbulence for the time being, in a magnetized disc 
with a purely toroidal magnetic field the magnetic tension tries 
to prevent radial motion associated with horseshoe U-turns im- 



posed by the planet's gravity, and attempts to force fluid ele- 
ment trajectories to follow field lines. When the magnetic field 
is strong enough, it is conceivable that magnetic tension can pre- 
vent horseshoe U-turns, cancelling the corotation torque. This 
might explain why the torque distribution in the 2D MHD simu- 
lations of Fromang et al . (2005), where/? = 2, reveals no corota- 
tion torque, but instead a magnetic torque with properties similar 
to those predicted by linear theory. However, a weak magnetic 
field, such as the one we consider here, cannot prevent horseshoe 
U-turns, as we have demonstrated in Sect. 15.11 In this case, we 
speculate that magnetic tension will lead to asymmetric horse- 
shoe U-turns, with downstream streamlines being closer to the 
corotation radius than the upstream streamlines. This is equiv- 
alent to saying that fluid elements in weakly magnetized discs 
will experience smaller changes to their angular momenta dur- 
ing horseshoe U-turns than they will in non-magnetized discs, 
with a consequence being that the density near the downstream 
separatrices will be reduced. Magnetic tension will thus reduce 
the positive corotation torque on the planet arising from fluid ele- 
ments undergoing inward horseshoe U-turns, and will reduce the 
magnitude of the negative corotation torque arising from fluid 
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elements undergoing outward horseshoe U-turns. This predicted 
behaviour is in broad agreement with the resuhs presented in 
Fig- US provided we ignore the fact that the modification to the 
corotation torque observed in the simulations is not quite sym- 
metric about the corotation radius. The origin of this offset about 
the corotation radius is not understood at present. A more de- 
tailed investigation of horseshoe dynamics in the presence of a 
weak toroidal magnetic field and MHD turbulence go beyond 
the scope of this work, however, and will be presented in a forth- 
coming publication. 

We finally comment that the results of this study have been 
obtained with rather large values of the disc aspect ratio [hp - 
0.1) and of the planet mass (Mp = 3 x lO^^^M*). Alth ough t his 
planet mass is relevant to type 1 migration (see Sect. I2.2.2l i. it 
is a priori not straightforward to adapt our results to smaller, 
more canonical values of hp and Mp, or to diff'erent values of the 
alpha viscosity parameter associated with the turbulent stress. 
This requires knowledge of the properties of the total corotation 
torque in laminar magnetized discs, and whether or not they are 
altered by MHD turbulence. 



5.3. Mass-adding procedure 



The MHD simulations presented in this study were carried 
out using the mass-adding procedure described in Sect. 12.2.11 
This procedure, which restores the initial density profile on a 
timescale set to 20 orbital periods, aims at obtaining a steady- 
state density profile, mass being otherwise continuously relo- 
cated by the turbulence towards the radial boundaries of the com- 
putational domain. Not including this procedure would make 
both the density profile and the running time-averaged tidal 
torque non-stationary quantities. In the simulations presented in 
Sect, m the choice for a restoring timescale of 20 orbital periods 
helps obtain a time-averaged density profile close to the initial 
profile, and which can be approximated as a power-law func- 
tion of radius over a rather large radial extent about the planet's 
orbital radius (see Fig. |7]i. This has enabled us to compare the 
torques of our simulations with analytical torque formulae de- 
rived in disc models with power-law density profiles. A different 
restoring timescale would result in a different density profile in 
a steady-state, and therefore in a diff'erent tidal torque. Although 
not presented here, we performed another MHD simulation with 
the parameters of Model 1 and a restoring timescale increased 
to 40 orbital periods. This additional run was performed with 
the NIRVANA code. We also carried out an additional 2D vis- 
cous simulation using the FARGO code, where the initial pro- 
files of the density and of the viscous alpha parameter corre- 
spond to their time-averaged counterpart in the MHD simula- 
tion (this is the same approach as that adopted in Sect. |4|i. The 
agreement between the laminar and turbulent simulations is very 
similar to that obtained with a restoring timescale of 20 orbital 
periods: the time-averaged turbulent and laminar torques agree 
to within ~ 10%, and the torque density distributions show the 
same features as inthe top-right panel of Fig. [T3] This agreement 
makes us confident that the mass-adding procedure does not in- 
troduce any other artifacts on the time-averaged tidal torque than 
the one naturally arising from changing the density profile in a 
steady-state by choosing a diff'erent restoring timescale. We fi- 
nally stress that this procedure should not be used in turbulent 
disc models involving a gap-opening planet. 



6. Conclusions 

In this paper we have investigated the tidal torque between a 
planet on a fixed circular orbit and its nascent protoplanetary 
disc, assuming that the latter sustains turbulence driven by the 
non-linear development of the MRl. Our main aim has been 
to investigate the properties of the corotation torque in such 
turbulent discs. In this initial study, we adopted simple 3D 
magnetized disc models with a locally isothermal equation of 
state, and neglected vertical stratification and non-ideal MHD ef- 
fects. We performed simulations using two different MHD codes 
(NIRVANA and RAMSES) over time spans of several hundred 
planet orbits. 

In one set of simulations, we considered the torque experi- 
enced by a planet orbiting at the outer edge of an inner disc cav- 
ity, in a region with a strong, positive density gradient. The run- 
ning time-averaged torque in this case was found to be strongly 
positive over several hundred planet orbits, in good agreement 
with the results of laminar viscous disc simulations whose vis- 
cous stresses were similar to that generated by the MHD tur- 
bulence. This positive torque strongly suggests the existence of 
an unsaturated corotation torque in weakly magnetized turbulent 
disc models. It implies that the presence of an inner cavity can 
stall type 1 migration, such that a planet trap can operate in tur- 
bulent discs. 

We also performed simulations of two disc models with 
power-law initial density profiles. We found that the running 
time-averaged torque of the turbulent disc models reached a sta- 
tionary value after 100 to 200 planet orbits, with decent agree- 
ment between it and a laminar viscous disc model with simi- 
lar viscous alpha parameter at the planet location. The relative 
difference between the averaged torques in turbulent and lami- 
nar disc simulations ranges from 20% to 60%, depending on the 
code and the disc model. We find remarkably good overall agree- 
ment between the results of the two codes, and given the scales 
of interest (the half-width of the planet's horseshoe region is a 
fraction of the pressure scale height), this agreement highlights 
the robustness of our results. 

Analysis of the velocity fields in our simulations shows 
the existence of horseshoe dynamics in turbulent discs with 
a weak toroidal magnetic field. The shape and width of the 
time-averaged horseshoe region much resemble those of lami- 
nar viscous discs with similar viscosity at the planet's location. 
Furthermore, close inspection of the averaged torque density 
distributions demonstrates that the differential Lindblad torque 
takes very similar values in turbulent and laminar viscous disc 
models, and there exists an unsaturated corotation torque in 
weakly magnetized turbulent discs. This analysis, however, also 
indicates the existence of an additional torque with respect to the 
non-magnetic case, which originates from the planet's horseshoe 
region. We argue that this additional torque is likely to be an ad- 
ditional corotation torque, possibly arising from the effects of 
magnetic tension. 

Although we have demonstrated conclusively the existence 
of unsaturated corotation torques in weakly magnetized turbu- 
lent discs, our study contains a number of simplifications that 
need to be addressed. We have considered disc and planet param- 
eters such that the horseshoe U-turn drift rate exceeds the ampli- 
tude of the turbulent velocity fluctuations. In the opposite case, 
the existence of horseshoe dynamics and a corotation torque is 
unknown. We have neglected vertical stratification and non-ideal 
MHD effects such as Ohmic resistivity, and it is fully expected 
that inclusion of these effects will modify the results presented 
here. We have also kept the planet on a fixed circular orbit and 
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simply measured the tidal torque as a time series. In reality the 
planet will migrate and experience evolution in its eccentricity 
and inclination as a result of interaction with the turbulent disc. 
We will address these issues in a series of forthcoming publi- 
cations, along with an in-depth investigation of the additional 
corotation torque in weakly magnetized discs that our simula- 
tions have uncovered. 
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